home *** CD-ROM | disk | FTP | other *** search
/ Libris Britannia 4 / science library(b).zip / science library(b) / MATH / MFLOAT10.ZIP / BESSEL.C < prev    next >
C/C++ Source or Header  |  1993-04-28  |  2KB  |  74 lines
  1. /*  This test program shows the advantages of mfloat numbers. */
  2. /*  The bessel function J0(x) is calculated by the series.    */
  3.  
  4. extern unsigned _stklen = 40000U;
  5.  
  6. #include "mfloat.h"
  7. #include <stdio.h>
  8. #include <string.h>
  9.  
  10. /**************************************************************************/
  11.  
  12. long double J0(long double x) {
  13.  
  14.   mfloat sum, sqrx, prod, mepsi;
  15.   long double epsi;
  16.   int i;
  17.  
  18.   epsi = 1E-20;
  19.   ldtomf(mepsi,epsi);
  20.   sqrm(ldexpm(ldtomf(sqrx,x),-1));
  21.   getonem(sum);
  22.   getonem(prod);
  23.   i = 0;
  24.   do {
  25.     i++;
  26.     divi(divi(multm(prod,sqrx),i),i);
  27.     if (i & 1) subm(sum,prod);
  28.     else addm(sum,prod);
  29.   } while (!gtm(mepsi,prod));
  30.   return(mftold(sum));
  31. }
  32.  
  33. /**************************************************************************/
  34.  
  35. long double J0x(long double x) {
  36.  
  37.   long double sum, sqrx, prod, epsi;
  38.   int  i;
  39.  
  40.   epsi = 1E-20;
  41.   sqrx = x*x/4;
  42.   sum = 1;
  43.   prod = 1;
  44.   i = 0;
  45.   do {
  46.     i++;
  47.     prod *= sqrx / i / i;
  48.     if (i & 1) sum -= prod;
  49.     else sum += prod;
  50.   } while (prod >= epsi);
  51.   return(sum);
  52. }
  53.  
  54. /**************************************************************************/
  55.  
  56. int main (void) {
  57.  
  58.   int accuracy;
  59.   long double x;
  60.  
  61.   printf("              Calculation of the bessel function J0(x)\n");
  62.   printf("              ========================================\n\n");
  63.   x = 100.0;
  64.   printf("mantissa words        result     (x = %5.1Lf)\n\n", x);
  65.   for (accuracy = 1; accuracy <= 15; accuracy++) {
  66.     setmantissawords(accuracy);
  67.     printf("%3i                J0(%5.1Lf) = %25.18Le\n", accuracy, x, J0(x));
  68.   };
  69.   printf("\n");
  70.   printf("IEEE arithmetic    J0(%5.1Lf) = %25.18Le\n\n", x, J0x(x));
  71.   printf("The accuracy of 12 mantissa words is needed to get a good result!\n");
  72.   return(0);
  73. }
  74.